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This  report  descrUies  the  derivation  of  a  set  of  ordinary  differential  equations  to 
model  radical  chain  polymerisation.  These  equations  have  the  mathematical 
property  of  stiffness  and  are  difficult  to  solve  numerically.  We  show  how  these 
equations  can  be  solved  efficiently  using  either  the  Gear  or  Kaps-Rentrop  method. 
We  also  show  how  the  kinetic  scheme  can  be  expanded  to  allow  for  the  presence  of 
contaminant  scavenger  molecules,  and  we  appiy  these  schemes  to  model 
experimental  results  for  the  polymeri-  -  uon  ofN-vinyl-2-pyrrolidone  obtained  from 
dilatometry  measurements. 
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Numerical  Solution  of  Stiff 
Ordinary  Differential  Equations 
for  Polymerisation  Kinetics 


1.  Introduction 

Materials  Research  Laboratory  is  currently  engaged  on  a  program  of  research  and 
development  on  a  number  of  cast-cured  polymer  bonded  explosives  (PBXs)  as  part 
of  the  the  development  of  an  Australian  Insensitive  Munitions  c^atality.  These 
fillings  are  formulated  by  dispersing  solid  explosive  filler  in  a  liquid  pre-polymer, 
adding  curing  agent,  mixing,  then  casting  into  the  munition  and  curing  in  situ. 
Among  the  PBX  formulations  being  studied  is  one  which  consists  of  RDX  in  an 
acrylic  binder  which  is  formed  in  situ  by  the  copolymerisation  of  2-ethylhexyl- 
acrylate  (EHA),  dioctylmaleate  (DOM),  and  N-vinyl-2-pynDlidone  (NVP). 

As  part  of  a  study  of  the  curing  of  the  binder  in  this  PBX  formulation  an 
investigation  of  the  mathematical  modeling  of  the  polymerisation  reactions  has  also 
started.  The  aim  is  to  use  the  modeling  as  an  aid  to  understanding  the  curing 
process  and  to  have  some  predictive  capability  for  cure  rate  as  a  fimction  of  various 
parameters  such  as  temperature,  inhibitors,  initiators  and  monomer  ratios.  We  have 
begun  by  considering  the  polymerisation  of  NVP  in  isolation  from  other  binder 
monomers  and  RDX.  The  polymerisation  is  initiated  using  2,2'-azobis 
(2-methylpropionitrile)  (AIBN)  and  the  rate  of  the  reaction  is  followed  by 
dilatometry. 

The  mathematical  model  of  this  process  is  described  by  a  set  of  coupled  ordinary 
differential  equations  (DDEs)  which  have  the  mathematical  property  of  stifhiess, 
which  means  that  the  dependent  variables  can  change  on  two  or  more  very  different 
scales  of  the  independent  variable.  This  seriously  complicates  the  numerical 
solution  of  the  coupled  equations  when  simple  explicit  numerical  schemes  are 
employed  because  the  stability  analysis  shows  that  very  small  time  steps  must  be 


used  even  in  regions  where  the  solution  is  only  slowly  varying.  For  the  NVP 
system  this  implies  up  to  one  million  iterations  of  the  solution  procedure  to  follow 
the  system  for  just  one  second. 

Fortunately  practical  algorithms  for  the  numerical  solution  of  stiff  ODEs  have 
been  devised  and  implemented  in  various  software  packages.  A  typical  examine  is 
the  method  due  to  Gear  [1],  which  is  implemented  in  the  NAG  routine  D02EBF  and 
is  available  at  MRL  [2].  The  solution  of  these  equations  is  so  error  prone,  however, 
that  it  is  both  comforting  and  advisable  to  have  a  completely  independent  check  on 
their  accuracy.  To  provide  such  a  practical  and  independent  method  of  checking 
the  solutions  of  our  kinetic  schemes  we  have  implemented  an  efficient  and 
relatively  simple  algorithm  for  the  numerical  solution  of  stiff  sets  of  ODEs.  It  is 
based  on  the  Kaps-Rentrop  method  [3]  as  described  by  Press  and  Teukolsky  [4]  and 
uses  an  implicit  Runge-Kutta  method.  The  algorithm  requires  no  more  than  a  few 
hundred  lines  of  FORTRAN  coding  and  could  easily  be  converted  to  BASIC  and 
run  on  a  PC  if  so  desired. 

The  purpose  of  this  report  is  firstly  to  derive  the  equations  which  model  the 
radical  chain  polymerisation  process,  and  then  to  describe  several  numerical 
techniques  for  the  efficient  solution  of  those  equations.  In  the  next  section  we 
describe  the  standard  kinetic  scheme  for  radical  chain  polymerisation  of  monomers 
such  as  NVP  and  then  briefly  discuss  some  approximate  solutions  of  these 
equations.  In  Section  3  we  then  describe  several  numerical  methods  for  their 
solution.  We  have  obtained  solutions  using  an  explicit  4th  order  Runge-Kutta 
routine,  the  NAG  routine  D02EBF  based  on  the  Gear  method,  and  solutions  using 
the  Kaps-Rentrop  method.  Using  the  Kaps-Rentrop  method  we  found  excellent 
agreement  with  previous  results,  but  with  far  less  computational  effort.  To  follow 
the  model  out  to  8  000  seconds  for  example  requires  only  100  solution  steps, 
whereas  to  check  the  results  using  the  explicit  Runge-Kutta  code  would  have 
required  approximately  one  billion  steps.  In  Section  4  we  extend  the  kinetic 
scheme  described  in  Section  2  to  include  the  presence  of  radical  scavenger 
molecules  and  illustrate  the  process  of  cure  inhibition/retardation.  We  also 
compare  the  numerical  solutions  with  experimental  data  on  the  polymerisation  of 
NVP  contaminated  with  an  as  yet  unidentified  impurity. 


2.  Kinetic  Scheme  for  NVP  Polymerisation 
and  Approximate  Solutions 


The  polymerisation  of  NVP  occurs  by  the  radical  chain  mechanism.  This  scheme  is 
described  in  detail  in  standard  texts  [5, 6],  and  is  briefly  summarised  here.  Radical 
polymerisation  is  a  chain  reaction  which  requires  the  steps  of  initiation, 
propagation,  and  termination.  The  initiation  step  involves  two  reactions;  ilist  an 
initiator  I  dissociates  to  form  a  pair  of  radicals  R* ,  and  then  in  the  second  step  a 
radical  combines  with  a  monomer  molecule  M  to  produce  a  chain  initiating  species 
M,*.  These  reactions  are  written  as 


I 


2R 


(1) 


R*  +M - ►Mj* 


(2) 


and  occur  with  rate  constants  and  kj  respectively  ^  In  the  propagation  step  the 
chain  radical  formed  in  the  initiation  step  (Mj  * )  grows  by  the  addition  of  very  large 
numbers  of  monomer  molecules  according  to  the  general  scheme 


n  successive 

M,  •  +  M - -  Mj* - -  (3) 

steps 


The  rate  constant  for  propagation  is  denoted  by  k^  and  its  numerical  value  is 
independent  of  the  size  of  the  growing  radical  after  the  first  few  additions.  The  final 
step  is  the  termination  of  polymer  growth  and  this  can  occur  in  one  of  two  ways, 
either  combination  or  disproportionation.  Both  processes  involve  a  bimolecular 
reaction  of  the  radical  sites  at  the  ends  of  growing  polymer  molecules.  A 
combination  reaction  results  in  a  single  molecule,  i.e.  two  ends  have  been  joined  by 
a  chemical  botxl  producing  a  larger  chain.  Disproportionation  involves  one  radical 
centre  gaining  a  proton  resulting  in  one  chain  with  a  saturated  chain  termination  and 
the  other  with  a  double  bond  termination.  Termination  can  also  occur  by  a 
combination  of  coupling  and  disproportionatioa  The  termination  step  is  generally 
represented  by 


M  •  +  M„* 

n  in 


M. 


(4) 


where  the  particular  mode  of  termination  is  not  specified,  and  the  rate  constant  k^  is 
the  sum  of  the  rate  constants  for  each  individual  termination  process. 

Equations  (1)  through  (4)  can  be  described  by  the  following  set  of  coupled  rate 
equations  for  the  concentrations  of  initiator  [I],  radicals  [R*  ],  monomers  [M],  and 
growing  polymer  molecules  [M*  ]. 


d[ll 

—  =  -kj[I]  (5) 


1  The  tcheme  depicted  in  equation  (l)Mmie  only  for  ■lymmetrical  initiator  L  A  more  typical  icheine  ia  ai  foUowi 
1 - -  Rl*  +R2*. 

with  only  one  of  the  decompoiiliao  producti  hein(  the  reactive  ipedea. 


T 


d[R*] 

dt 


d[M] 

dt 


d[M*] 

dt 


=  2kj[I]-ki[R*][M] 


=  -kj[R*][M]-kp[M*][M] 


=  kj[R*][M]-2k,[M*^ 


(6) 


(7) 


(8) 


Equation  (S)  indicates  a  simple  exponential  decay  for  the  initiator  concentration  [I], 
and  so  the  set  of  equations  (S)  through  (8)  can  be  replaced  by  the  more  convenient 
set 


dR* 

dt 

=  2fkjlQexp(-kjt)-kjR*M 

(9) 

dM 

dt 

=  -kjR’M-kpMM* 

(10) 

dM* 

dt 

=  kjR*M-2kj(M*)^ 

(11) 

where  we  have  introduced  the  initiator  efficiency  f,  which  is  defined  as  the  fraction 
of  the  radicals  produced  in  the  decomposition  reaction  which  initiate  polymer 
chains.  The  value  of  f  is  usually  less  than  one  and  we  have  used  a  value  of  0.47, 
which  was  detennined  by  Braun  and  Quella  [7]  for  the  initiation  system  being  used. 
We  have  also  dropped  the  bracket  notation  and  simply  let  R* .  M  and  M*  denote  the 
concentrations  of  radical,  monomer,  and  growing  polymer  molecules  respectively. 

The  rate  constants  k^,  k.,  and  k,  will  determine  the  degree  of  stiffness  of  the 
set  of  equations  (9)  through  (1 1).  The  rate  of  decomposition  of  the  initiator  k^ 
(equation  1)  can  to  a  large  degree  be  selected  by  choosing  the  appropriate  chemical 
system,  and  the  rate  of  the  overall  polymerisatitm  process  can  also  be  adjusted 
substantially  by  adjusting  the  concentration  of  the  initiator  and  conditions  to  suit  the 
puipose  of  the  polymer  production  process.  For  example,  in  a  compilation  of 
decomposition  rate  constants  edited  by  Brundip  and  lmmergut[8],  k^  varies  over  a 
range  from  10"’°  s'*  to  10'*  s  *  for  various  initiators  under  various  conditions.  In 
the  case  of  the  PBX  the  situation  is  complicated  by  the  use  of  two  peroxide  initiators 
as  well  as  cobalt  (II)  bisacetylacetonoate  (CoAA).  The  CoAa  allows  the  PBX  to 
cure  at  room  temperature  because  it  facilitates  the  peroxide  dissociation,  i.e.  k^  is 
increased  by  the  presence  of  CoAA.  We  have  not  yet  measured  k^  for  the 
peroxides  used  in  the  PBX  in  the  presence  of  CoAA,  and  for  this  study  of  the  NVP 
polymerisation  a  simpler  initiation  system  using  AIBN  was  chosen  so  that 
comparisons  could  be  made  with  literature  data.  Braun  and  (^eUa[7]  have 


10 


measured  and  found  a  value  of  1.62  x  10'^  s'^  at  60’C  for  AIBN.  Our 
experiments  have  not  yet  yielded  a  value  for  k^. 

The  rate  constants  for  the  propagation  step  k^  (equation  3)  and  the  termination 
step  k,  (equation  4)  have  been  extensively  studied  for  large  numbers  of  monomers. 
Some  of  the  data  relevant  to  monomers  in  the  PBX  binder  have  been  found  in  the 
literature,  but  the  relevant  data  is  very  limited.  However  it  is  known  that  for  most 
systems  the  values  of  kp  and  k^  lie  in  the  range  10^  to  10^  1  mole'*  s'*  and  10®  to 
10*  1  mole'*  s'*  respectively  [6]. 

The  rate  constant  k.  for  the  initiation  step  (equation  2)  is  not  known  and  rate 
constants  for  reactions  of  this  type  have  not  been  found  in  the  literature.  However 
the  rates  of  reactions  involving  the  so  called  "primary  radicals"  and  the  monomer 
should  be  of  similar  magnitude  to  the  propagation  rate  constant  kp  (equation  3). 
This  lack  of  knowledge  regarding  the  value  of  k.  is  not  critical  as  we  will  shortly 
demonstrate  that  it  has  little  effect  on  the  overall  rate  of  polymerisation. 

With  the  above  considerations  in  mind  we  have  chosen  to  use  the  following  set 
of  values  for  the  four  rate  constants  for  our  preliminary  investigation  of  the  system: 

kj  =  1.6  X  10'^  s'* 
k.  =  1.0  X  10^  1  mole  *  s'* 
kp  =  1.0  X  10^  1  mole'*  s'* 
k,  =  1.0  xio’ I  mole  *  s'* 


The  range  of  values  for  the  rate  constants  noted  above  means  that  the  set  of 
equations  (9)  through  (1 1)  will  have  the  mathematical  property  of  stiffness. 
Technical  definitions  of  stiffness  can  be  given  in  terms  of  the  eigenvalues  of  the 
Jacobian  matrix  formed  from  the  equations  [1],  but  here  it  simply  suffices  to  note 
that  any  set  of  equations  in  which  the  dependent  variables  can  change  on  two  or 
more  very  different  scales  of  the  independent  variable  are  called  stiff.  Stiff  sets  of 
equations  are  notoriously  difficult  to  solve  numerically  using  explicit  schemes 
because  the  stability  of  the  scheme  is  governed  by  the  size  of  the  time  step  needed 
to  resolve  details  of  changes  occurring  on  the  fastest  time  scale.  This  means  that 
very  small  time  steps  are  required  to  follow  changes  over  very  much  longer  periods 
of  time,  even  though  the  solutions  may  not  be  varying  rapidly  on  this  time  scale. 
Before  describing  a  practical  method  for  overcoming  this  problem  in  the  next 
section,  we  first  briefly  describe  an  tqrproximate  solution  of  the  equations  (9) 
through  (1 1)  to  get  a  feel  for  their  behaviour. 

The  analysis  of  radical  chain  polymerisation  in  most  textbooks  begins  by 
assuming  the  steady  state  condition,  i.e.  the  assumption  is  made  that  the  rate  of 
initiation  is  equal  to  the  rate  of  termination,  so  that  the  concentration  of  free  radicals 
becomes  essentially  constant  very  early  in  the  reaction.  We  use  this  assumption 
ourselves  in  a  moment,  but  first  note  that  the  decrease  in  monomer  concentration 


during  a  typical  polymerisation  of  the  type  involved  in  our  PBX  binder  studies  is 
very  slow,  decreasing  by  approximately  10%  over  a  time  of  several  hours.  This 
means  that  we  can  replace  the  variable  M  in  equation  (9)  by  its  initial  value  at  t 
equals  zero.  Also,  the  value  of  is  such  that  the  exponential  in  equation  (9)  may 
be  replaced  by  unity  for  times  shorter  than  one  hour.  Equation  (9)  then  becomes 


dR* 

dt 

with  solution 


2fkJ  -k.M  R* 
do  10 


(12) 


R‘(t) 


— i^{l-exp(-kjM  t)} 


(13) 


indicating  a  rise  to  a  steady  state  value  of  2fkjljj/kjMjj  in  a  time  of  order  (k-M^)''. 
Appropriate  experimental  values  for  and  are  9.08  moles  r‘  and  1.52  x  10"^ 
moles  r'  respectively  and  these  lead  to  a  steady  state  value  of  R*  =  2.5  x  10'^^ 
moles  r*  in  a  time  of  about  10"^  s. 

We  now  set  dM*  /dt  equal  to  zero  in  equation  (11)  and  assume  that  the  radical 
concentration  has  reached  its  equilibrium  value  R*  .  We  then  obtain  the  following 
expression  for  the  equilibrium  value  of  M* 


eq  0 


(14) 


Using  the  values  for  the  constants  given  above  this  leads  to  a  value  of 
M*  =  3.4  X  10  *  mole  T*.  We  could  solve  equation  (11)  exactly  to  obtain  an 
estimate  of  the  time  taken  for  M*  to  reach  this  equilibrium  value  by  substituting  the 
solution  for  R*  given  by  equation  (13)  into  equation  (11).  The  resulting  equation 
is  a  Riccati  equation,  which  is  not  easy  to  solve,  but  solutions  can  be  expressed  in 
terms  of  the  solutions  of  two  auxiliary  equations  [9].  We  have  not  pursued  this 
tq)proach  however  as  we  know  from  prior  experience  that  this  time  scale  is  on  the 
order  of  seconds  or  less. 

We  now  consider  the  decrease  in  moromer  concentration  on  the  slower  time 
scale.  Assuming  that  both  R*  and  M*  do  not  vary  greatly  from  their  equilibrium 
values  as  M  begins  to  decrease,  we  can  treat  equation  (10)  as  a  simple  equation  in 
the  single  variable  M.  The  solution  is  then 

M(t)  =  Mj,cxp(-kt)  (15) 

whwe 


k  =  kjR*  +k  M* 

•  eq  P 


eq 


(16) 
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With  the  values  for  Iq.k^.R*  ^M*  used  above  k  can  be  approximated  very 
closely  by  k^M*  arid  the  ex^nmial  c^  be  expanded  to  first  order  to  give 

M(t)  =  Mj,  { 1  -  kp  [(kj/k,)  f  y  t}.  (17) 


indicating  that  a  plot  of  monomer  concentration  versus  time  for  the  early  stage  of 
monomer  consumption  should  lie  on  a  straight  line  and  have  a  negative  gradient  of 
Mokp  {(kj/kj)  f  Iq}  That  this  is  indeed  the  case  can  be  seen  in  Ftgure  1,  which 
shows  a  relationship  between  monomer  concentration  and  time,  derived  from 
experimental  data  obtained  with  a  dilatometer  during  an  NVP  polymerisation 
reaction.  The  monomer  concentration  versus  time  data  was  obtained  from  the 
volume  change  versus  time  data  by  using  a  conversion  factor  which  incorporates  the 
densities  of  NVP  and  the  polymer  at  the  temperature  of  the  reaction.  This 
conversion  factor  is  only  known  at  20*C  at  present  and  this  value  was  used  even 
though  the  data  was  obtained  at  60‘C.  The  data  cannot  therefore  yet  be  used  to 
obtain  accurate  values  of  kinetic  parameters,  but  it  can  be  used  to  provide  a 
comparison  for  the  model  output  and  approximate  values  of  parameters. 
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Figure  I :  Experimental  data  for  monomer  concentration  versus  time  for  NVP. 


Given  that  k^,  f,  1  and  are  known,  a  value  for  the  ratio  kp^/k,  can  be  obtained 
from  the  straight  line  fit  to  the  experimental  data  in  Figure  1 .  In  practise,  a  more 
accurate  method  ui  doing  this  is  to  obtain  the  gradient  dM/dt  for  several  values  of 
the  initial  initiator  concentration  y  From  equation  (17)  we  easily  have 


-dM/dt  =  [Mj,(fk^>«(kp/k,)>«]Ij,>« 


(18) 


and  so  a  plot  of  gradient  versus  again  leads  to  a  value  of  Measurements 
of  as  a  ftmction  of  temperature  lead  to  useful  information  about  die 
thermochemistry  of  polymerisation,  but  it  is  importani  to  note  that  values  for  and 
k,  separately  can  only  be  obtained  by  non  steady-state  methods  [S] .  We  ai^ 
currently  reviewing  the  need  to  use  a  non-steady-state  method  such  as  the  rotating 
sector  method  [6]  to  obtain  the  additional  experimental  data  required  to  separate  k^ 
and  k^. 

Equation  (17)  also  indicates  that,  with  the  range  of  values  for  the  rate  constants 
we  are  using,  the  initial  rate  of  polymerisation  is  effectively  independent  of  the  rate 
constant  k.  and  justifies  our  assertion  that  knowledge  of  an  exact  value  for  k.  is  not 
required  for  the  polymerisation  experiments  we  are  planning. 


3.  Numerical  Solution  of  Kinetic  Equations 


The  initial  rapid  rise  of  R*  to  its  equilibrium  value  R*  on  a  time  scale  of  the 
order  of  milliseconds  is  easily  demonstrated  using  a  sim^e  4th  order  Runge-Kutta 
integration  routine.  Figures  2  and  3  show  solutions  of  equations  (9)  through  (11) 
obtained  using  the  program  RKUTTA.  (The  source  code  for  the  programs 
referred  to  in  this  report  arc  available  from  the  authors).  The  rate  constants  kj,  kj, 
kp,  kj  and  the  initial  constants  and  1^  were  set  to  the  values  used  in  the  previous 
section.  Figure  2  was  calculated  using  a  time  step  At  of  1  ps  and  shows  that  R* 
reaches  an  equilibrium  value  of  2.5  x  10'*^  mole  T*  in  a  time  of  approximately 
0.5  ms,  which  is  in  good  agreement  with  the  approximate  calculation  in  the  previous 
section.  On  this  time  scale  M  remains  constant  while  M*  increases  linearly  with 
time,  as  can  be  seen  in  Figure  3,  which  shows  both  R*  and  M*  between  t  =  0  and 
5.0  ms  calculated  using  a  time  step  of  5.0  ps. 

To  use  the  program  RKUTTA  to  follow  the  time  variation  of  either  M*  or  M 
would  be  completely  impractical.  A  time  step  of  the  order  of  10  ps  would  be 
needed  to  establish  the  equilibrium  value  of  R*  accurately,  and  then  many  millions 
of  steps  would  be  needed  to  track  either  M*  or  M.  The  sensible  procedure  to 
follow  for  the  solution  of  stiff  sets  of  equations  such  as  (9)  through  (1 1)  is  to  use 
one  of  the  software  packages  specifically  designed  for  the  solution  of  these 
equations.  These  are  usually  based  on  implicit  schemes  with  variable  step  length 
and  automatic  local  error  control.  A  particularly  well  known  procedure  for  the 
solution  of  stiff  sets  of  equations  is  the  Gear  method  [1],  which  is  available  at  MRL 
as  subroutine  D02EBF  of  the  NAG  group  of  software  packages.  We  have  used  the 
Gear  method  to  follow  the  full  time  dependence  of  the  variables  in  equations  (9) 
through  (11). 
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Figure  2 :  Radical  concentration  versus  time  calculated  from  program  RKUTT 4 
with  Ar  =  /  jis. 


Time,  ms 


Figure  3:  Concentrations  cf  radicals  and  growing  polymer  molecules  calculated 
from  program  RKUTTA  with  At  =  5.0  MJ. 
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Figure  4  shows  the  time  dependence  of  M*  on  a  time  scale  of  seconds.  We  see  that 
M*  approaches  a  steady  state  value  of  3.4  x  10  *  moles  1*,  which  is  again  in  good 
agreement  with  the  approximate  solution  in  the  previous  section.  Figure  5  shows 
R* ,  M*  and  M  over  a  longer  time  scale.  We  see  that  both  R*  and  M*  vary  little 
from  their  equilibrium  values,  justifying  the  assumption  made  in  the  previous 
section,  and  that  over  this  time  scale  M  decreases  linearly  with  time,  again  in 
agreement  with  the  analysis  in  the  previous  section.  To  display  the  exponential  time 
dependence  of  M  requires  following  the  solution  over  a  much  longer  time  scale. 

We  have  done  this  and  show  the  solutions  in  Figure  6.  The  exponential  decay  of  M 
is  now  clearly  evident. 


Tine.  3 


Figure  4:  Concentration  of  growing  polymer  molecules  versus  time  calculated 
using  the  Gear  method. 


We  have  used  the  subroutine  D02EBF  to  check  some  of  the  assumptions  and 
conclusions  of  the  previous  section.  In  particular,  our  analysis  predicts  that  over  a 
time  scale  of  just  a  few  hours  M  should  decrease  linearly  with  time,  and  that  the 
gradient  of  this  time  should  be  independent  of  k.,  and  depend  on  k^  and  k^  only  in 
the  ratio  kp^/k,.  We  have  checked  these  results  by  varying  kj  over  the  range  1x10^ 
to  1  X  lO'*  and  found  that  it  has  had  no  effect  on  the  time  dependence  of  M  over  the 
time  scale  of  interest.  We  have  also  varied  the  values  of  kp  and  k^,  and  again  found 
that  the  time  dependence  of  M  is  unaltered  if  the  value  of  the  ratio  kp^/kj  remains 
constant. 
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Figure  5:  Concentration  of  radicals,  monomer,  and  growing  polymer  chains 
versus  time  calculated  using  the  Gear  method. 


Figure  6:  Concentration  of  radicals,  monomer,  and  growing  polymer  chains 
versus  time  calculated  using  the  Gear  method. 
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The  Gear  method,  as  implemented  in  IX^EBF,  is  both  simple  and  efficient  to  use, 
and  capable  of  providing  solutions  to  any  desired  degree  of  accuracy  by  setting  an 
appropriate  value  of  the  parameter  TOL  in  the  calling  sequence.  It  is  however  a 
"black  box”,  and  on  the  few  occasions  when  die  method  has  failed  it  has  not  been 
easy  to  trace  the  underlying  cause  of  this  failure.  To  overcome  this  problem,  and 
also  to  provide  a  completely  independent  means  of  checking  the  accuracy  of  the 
solutions  produced  by  D02EBF,  we  have  recently  implemented  a  new  method  for 
the  solution  of  stiff  sets  of  equations.  Hie  method  is  due  to  Kaps  and  Rentrop  [3], 
and  the  implementation  described  here  was  devised  by  Press  and  Teukolsky  [4]. 

A  set  of  non-linear  ODEs  can  be  written  in  the  form 


y  =  f(y) 


(19) 


where  y  is  the  vector  (y^,  y^, . y^)  of  N  dependent  variables,  and  the  prime 

denotes  differentiation  with  respect  to  time.  The  Kaps-Rentrop  method  seeks  a 
solution  of  the  form 

yCto+h)  =  yo+  SC.kj  (20) 

i=l 


where  the  corrections  are  found  by  solving  s  linear  equations  of  the  form 


(l-ThfO-k.  =  hf(y^+Za^kp 


+  hr* 


i-1 

Z  Yii  k: ,  I  =  1, 


s. 


(21) 


Here  the  coefficients  y,  Cj,  and  Yy  are  fixed  constants  independent  of  the 
problem.  Automatic  step  size  adjustment  is  provided  using  the 
Runge-Kutta-Fehlberg  method  [4];  two  estimates  having  the  form  of  equation  (20) 
are  computed,  one  of  higher  order  than  the  other;  the  difference  between  the  two 
then  leads  to  an  estimate  of  the  local  truncation  error,  which  can  then  be  used  for 
step  size  control.  Control  of  the  local  step  size  error  can  be  maintained  by  a 
suitable  choice  of  the  parameter  EPS  in  the  calling  program. 

Tables  1  through  3  illustrate  the  degree  of  accuracy  obtainable  using  both  the 
Gear  and  Kaps-Rentrop  methods.  It  should  be  noted  that  the  Gear  method  is 
implemented  in  double  precision,  while  the  Kaps-Rentrop  method  uses  only  single 
precision.  The  KAPS  program  is  marginally  faster  than  the  GEAR  program,  but 
both  take  no  more  titan  a  few  seconds  of  CPU  time. 
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Table  1 :  Gear  Method :  TOL  =  W* 


t(s) 

R*  X  10'2 

M 

• 

X 

o 

00 

500 

2.5364 

8.9412 

3.3680 

1000 

2.5587 

8.7939 

3.3540 

2000 

2.6032 

8.5048 

3.3789 

3000 

2.6484 

8.2274 

3.3007 

4000 

2.6935 

7.9612 

3.2744 

Table  2  :  Gear  Method :  TOL  =  70  '“ 


t(s)  R*  X  10'2  M  M*  X 10* 


500 

2.5401 

8.9284 

3.3674 

1000 

2.5625 

8.7796 

3.3540 

2000 

2.6075 

8.4911 

3.3272 

3000 

2.6526 

8.2143 

3.3007 

4000 

2.6978 

7.9487 

3.2744 

3:  Kaps-Rentrop  Method :  EPS  -  10* 

t(s)  R*  X 10'^  M 

M*  X  10* 

500 

2.5432 

8.9284 

3.3689 

1000 

2.5676 

8.7796 

3.3568 

2000 

2.6113 

8.4911 

3.3291 

3000 

2.6578 

8.2143 

3.3035 

4000 

2.7021 

7.9457 

3.2766 
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4,  Extended  Kinetic  Scheme  and  Numerical  Solution 


The  analysis  presented  in  Section  2  for  radical  chain  polymerisation  using  typical 
values  for  the  rate  constants  ^  and  the  initial  constants  and  cleaily 

shows  that  the  monomer  concentration  will  decay  lineariy  with  time  for  times  less 
than  one  or  two  hours.  Slight  variations  in  some  of  the  values  of  these  constants 
will  not  change  this  general  picture,  the  slope  of  the  line  may  change,  but  the 
decrease  will  still  be  linear  with  time.  Figure  7  shows  some  experimental  points  for 
NVP  monomer  concentration  versus  time  which  were  obtained  dilatometrically  in 
the  same  manner  as  described  in  Section  2.  The  only  significant  difference  in  the 
experimental  procedure  was  that  a  different  batch  of  NVP  was  used.  The 
pronounced  downward  concavity  of  these  points  indicates  that  some  mechanism 
other  than  simple  radical  chain  polymerisation  must  be  operating. 

To  explain  these  results  we  have  assumed  the  existence  of  a  contaminant 
scavenger  molecule  which  removes  both  radicals  and  activated  monomer.  We  can 
include  the  presence  of  such  a  contaminant  by  slightly  expanding  the  kinetic  scheme 
outlined  in  section  2.  We  define  a  scavenger  concentration  [S  *  ]  which  acts  on  both 
[R*  ]  and  [M*  ]  as  follows 


R*  +S*  — 

—  RS 

(22) 

M*  +S*  — 

—  MS 

(23) 

Formation  of  the  molecular  species  RS  and  MS  removes  both  R*  and  M*  from 
active  participation  in  the  radical  chain  growth  process.  We  assume  rate  constants 
kj  and  k^  for  equations  (22)  and  (23)  respectively.  Equations  (22)  and  (23)  can 
now  be  included  in  an  expanded  kinetic  scheme  which  has  the  following  form 

d[I] 

—  =  -kj(I]  (24) 

dt 

d[R*] 

-  =  2k^[I]-kj[R*][M]-k,[R’]IS*]  (25) 

dt 


d[M] 

-  =  -kiIR*l[M]-k  [M*KM] 

dt  ‘  P 


d[M*  ] 

- -  *  ki[R*][M]-2k,[M*]2-k2(M*nS*] 


dtS*] 


(26) 


(27) 


dt 


-k,  [R*][S*]-kj[M*][S*] 


(28) 


We  again  simplify  this  set  of  equations  by  incoiporating  the  exponential  solution  of 
equation  (24)  and  dropping  the  bracket  notation.  The  equations  then  become 


dR* 

dt 

=  2fkjI^jexp(-kjt)-kjR*M-kjR*S* 

(29) 

1  dM 

1  — 

' 

=  -kjR*M-kpM*M 

(30) 

'  dM* 

dt 

=  kjR*M-2k,(M*)2-k2M*S* 

(31) 

dS* 

i 

dt 

=  -k,R*S*  -k2(M*S*) 

(32) 

i 


We  have  solved  the  set  of  equations  (29)  through  (32)  using  both  the  NAG  routine 
D02EBF  and  the  Kaps-Rentrop  method  and  found  the  results  to  be  identical  within 
the  error  tolerances  of  both  schemes.  We  have  no  knowledge  of  the  values  of  kj 
and  k^  so  we  have  simply  set  kj  =  kj  =  k,  and  used  k  as  a  parameter  to  be  varied 
until  agreement  can  be  found  with  the  experimental  results.  The  curve  shown  in 
Figure  7  was  calculated  using  k  =  5.0  x  10^,  with  the  other  constants  having  their 
previous  values.  The  good  agreement  with  experiment  indicates  that  our 
assumption  of  the  presence  of  a  contaminant  scavenger  molecule  appears  to  be 
correct. 


5.  Conclusion  and  Discussion 

This  report  has  described  our  analysis  of  the  kinetics  of  the  polymerisation  of  NVP. 
We  have  discussed  the  kinetic  scheme  appropriate  to  radical  chain  polymerisation 
and  derived  appropriate  equations  to  model  this  system.  We  have  shown  how  to 
obtain  approximate  solutions  to  these  equations,  and  also  discussed  appropriate 
numerical  algorithms  for  the  efficient  solution  of  these  equations.  Our  analysis  has 
shown  that  by  following  the  polymerisation  rate  using  dilatometry  methods  we  can 
obtain  a  value  for  the  ratio  k^^/k,,  but  that  further  experimental  information  is 
required  to  obtain  values  for  k^  and  k,  separately.  We  are  currently  evaluating  the 
feasibility  of  several  experimental  techniques  to  provide  this  information. 

We  have  applied  the  mathematical  model  to  some  preliminary  experiments  on 
NVP  polymerisation  and  found  good  agreement  between  the  model  and  the 
experimental  results.  We  plan  to  expand  the  experimental  part  of  this  program 
considerably,  the  objective  being  to  determine  the  kinetic  parameters  necessary  to 
model  a  terpolymerisation  involving  the  three  monomers  NVP.  EHA  and  DOM.  In 
addition  we  wish  to  study  the  effects  of  inhibitors  (radical  scavanger  molecules). 
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Figure  7:  Comparison  of  experimental  values  of  monomer  concentration  versus 
time  with  theoretical  prediction  from  extended  polymerization  scheme. 

We  anticipate  that  the  model  we  have  derived  for  NVP  polymerisation  wiU  be 
applicable,  with  slight  modifications,  to  the  additional  chemical  systems  we  plan  to 
investigate,  and  expect  that  both  the  Gear  and  Kaps-Rentrop  algorithms  will  provide 
efficient  and  accurate  solution  ptocedures  for  the  analysis  of  future  experimental 
results. 
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